JLC 


AD-A204  404 


SHlS) 


4*.  '.AMS  -zr  »e"»Of'MlNG  OWGANIiATIQN 

University  of  Minnesota 


Ita.  Office  svmsol 

fif  ap0<itfad<#> 


ACOMCSS  iCity.  SIM  ana  X.IP  C-ioai 

1919  University  Avenue 
St.  Paul,  MN  55104 


Approved  for  public  release; 
distribution  unlimited. 


5.  MONITORING  ORGANIZATION  RSRORT  NuMSCRiSI 

AFOSR  TR.  8  9-  001  5 


7*.  NAMC  OP  MONITORING  ORGANIZATION 

Air  Force  Office  of  Scientific  Research/NM 


76.  AOORCSS  iCity.  Sian  ana  ZIP  Coaai 

Building  410 

Bolling  AFB,  DC  20332-6448 


NAME  0^  IIUNQING.SRONSORiNG 
ORGANIZATION 

AFOSR/NM 


.  OFf  icfi  svMaow 

fif  oppiicoduj 

NM 


9.  RROCURSMENT  instrument  identification  number 

AFOSR-87-0268 


8c.  aOCRESS  idly.  Staia  ana  ZIP  Coaai 

Building  410 

Bolling  AFB  DC  20332-6448 


10.  SOURCE  OF  FUNDING  NOS. 


I?  ftTLs 

"Local  and  Global  iGiJhn'tjotii  Fo.'" 


frqcram 
element  no. 


61102F 


froject 

TAS* 

UNtT 

NO. 

NO. 

^0 

2304 

A9 

!I  -mc|CU.t<-\Cy  oV  Per\OC\\C  ^Q^O+-|or\s  04  f-c,  rnv(2.-te  r -■  CepN^rvL-VAA 

y  Fonc-VvDnCyV  cTctoCt-horyj 


12.  FERSONAL  AUTnORISI 

Stech.  Harlan 


1  13a  tyfe  of  refort 

136.  Time  covered 

14.  OATE  OF  REFORT  /Yr.  Ma..  Oayi 

IS.  face  count  i 

1  FetiaI 

FROM  870715  Ta880915 

881118 

—I— _ 

18.  SUFFLEMENTaR  V  notation 


17  ccsati  cooes 

t8.  SU64CCT  TCP^MS  tCamtmu^  on  pvi/crat  if  nactumry  ana  lOonttfy  fry  Moca  aymbtri 

F'ElO  '  GROUF  SUB.  OR 

Numerical  Methods,  Functional  Differential  Equations, 

Hopf  bifurcation.  Periodic  Solutions 

1 

' 

1  i9.  aaSTRACT  Cjfttinuw  an  naetumry  and  lOtfin/y  9y  04o€a  numoarp  1 

^^This  project  was  a  continuation  of  an  ongoing  study  of  numerical/analytic  techniques 
for  the  identification  of  periodic  solutions  to  functional  differential  equations. 
The  techniques  developed  apply  to  very  general  classes  of  equations,  and  have  been 
Implemented  on  a  variety  of  specific  model  problems. 


The  areas  investigated  involve  techniques  and  information  not  attainable  by  standard 
simulation  methods.  The  work  completed  can  roughly  be  subdivided  according  to  the 
local  (Hopf  bifurcation)  analysis  in  the  neighborhood  of  equilibria,  and  global  tracking 
methods  for  following  1-parameter  families  of  periodic  orbits  and  examining  their 


secondary  bifurcation  structure. 


Za.  OlSTRlSUTICN/AvAlLASILlTV  OF  ABSTRACT 
UNC’.ASSIFiEO/UNLiMiTEO  S  SAMBAS  RFT.  "Z.  OTIC  USERS  C 


22m.  name  of  RESFCnSiBlS  iNOiviOual 


21  abstract  security  classification 

I  MA  C  Vc 


77^  rCuIF^CNf  NVJMtCA 
'faciyA* 

(202)  767-5028 


ill 


a 


22e.  OFFICE  SYMBOL 

NM 


edition  of  1  JAN  73  IS  obsolete. 

19 


lUHlinittlElt 


S  *AGc 


Dr^Air^^^achman^ 
00  FORM  1473,  33  APR 


SECL 


AhOSR  TK.  69-0015 


Final  Technical  Report 
on 


LOCAL  AND  GLOBAL  TECHNIQUES  FOR  THE  TRACKING 
OF  PERIODIC  SOLUTIONS  OF  PARAMETER-DEPENDENT 
FUNCTIONAL  DIFFERENTIAL  EQUATIONS 

AFOSR-87-0268 

July  15,  1987  -  September  15,  1988 

(AFSC) 

■i-'  £nc!  is 
••x.,  Ai-r^  190-12. 

-ri  informstior!  Division 


Accession  FoT 

HTIS  GRAStl 
OTIC  tab 

Unanno'inood 

Justi 


Distribution/ 

Avallablllt 

y  Codas 

md/er 

Lai 

• 

■ 

Dlst 

aJ 

Avail  t 
Speci 

I 


Harlan  w.  Stech 
Department  of  Mathematics 
and  Statistics 

University  of  Minnesota  -  Duluth 
Duluth,  Minnesota  55812 

18  November,  1988 


< 


Prepared  for 


Air  Force  Office  of  Scientific  Research 
Dr.  Arje  Nachman,  Program  Manager 
Building  410 

Bolling  AFB  DC  20330-6448 


B9  2  lO  112 


SUMMARY 


This  project  was  a  continuation  of  an  ongoing  study  of 
numerical/analytic  techniques  for  the  identification  of  periodic 
solutions  to  functional  differential  equations.  The  techniques 
developed  apply  to  very  general  classes  of  equations,  and  have  been 
implemented  on  a  variety  of  specific  model  problems, 

"Local"  techniques  refer  to  methods  that  apply  to  the  problem  of 
analyzing  the  Hopf  bifurcation  structure  of  small  periodic  orbits 
of  multiparameter  systems.  A  FORTRAN  code,  BIFDE,  was  modified  to 
allow  analysis  of  certain  nongenetic  bifurcations  in  general 
systems  with  infinite  delay.  Related  questions  have  been 
invest igated . 

"Global"  tracking  methods  have  been  developed  to  study  the  growth 
and  parameter  dependence  of  global  Hopf  bifurcations . 

Investigations  have  centered  on  the  development  of  spline-based 
approximation  techniques  and  their  implementation  in  a  FORTRAN  code 
FDETRAK.  The  capabilities  of  the  code  have  been  expanded,  and 
applied  to  scalar  and  two-dimensional  delay  difference  systems. 


RESEARCH  OBJECTIVES 


This  project,  along  with  the  grant  AFOSR-86-0071  with  the  same  name 
(budgeted  through  the  Virginia  Polytechnic  Institute)  initiated  the 
a  thorough  study  of  the  use  of  numerical  techniques  in  the  analysis 
of  general  parameter -dependent  functional  differential  systems. 

The  research  of  this  grant  was  simultaneously  supported  by  the 
National  Science  Foundation  grant  DMS-8701456,  and  a  grant  from  the 
DARPA  Program  in  Applied  and  Compuatational  Mathematics 
(subcontracted  through  Michigan  State  University.) 

The  areas  investigated  involve  techniques  and  information  not 
attainable  by  standard  simulation  methods.  The  work  completed  can 
roughly  be  subdivided  according  to  the  local  (Hopf  bifurcation) 
analysis  in  the  neighborhood  of  equilibria,  and  global  tracking 
methods  for  following  1-parameter  families  of  periodic  orbits  and 
examining  their  secondary  bifurcation  structure. 

At  this  early  stage  of  the  research,  emphasis  was  placed  on  the 
establishment  of  the  numerical  characteristics  of  some  of  the 
proposed  methods.  An  algorithm  of  the  PI  for  the  analysis  of  the 
bifurcation  structure  of  Hopf  bifurcations  was  chosen  to  be 
implemented  in  a  general  purpose  code  BIFDE.  A  graduate  student, 
Archana  Sathaye  was  supported  by  AFOSR-8 6-0071  to  assist  in  the 
implementation.  Her  work  was  continued  by  a  graduate  student, 

Nihal  Aboud,  under  the  support  of  this  grant.  Alternate  numerical 
implementations  (symbolic  manipulation)  and  techniques  for 
obtaining  the  necessary  bifurcation  data  were  addressed,  as  well. 

Complimentary  to  an  understanding  of  the  local  (Hopf  bifurcation) 
structure  of  multi-parameter  FDE  is  the  use  of  global  tracking 
methods  to  gain  insight  into  the  secondary  bifurcation  structure  of 
these  problems.  A  goal  of  this  study  has  been  to  design  a  general 
purpose  FORTRAN-based  code  for  the  analysis  of  scalar 
delay-difference  equations,  as  well  as  certain  multi-dimensional 
generalizations.  The  code  has  been  tested  on  a  wide  variety  of 
scalar  equations  and  a  system  modeling  "chugging”  in 
liquid-propellant  fuel  rockets. 


RESEARCH  STATUS 


Concerning  the  local  problem,  various  extensions  were  made  on  an 
Liapunov-Schmidt  based  algorithm  of  the  PI.  The  algorithm,  which 
leads  to  a  direct  means  of  determining  the  stability  and  direction 
of  bifurcation  in  functional  differential  equations  requires 
substantial  computation.  A  FORTRAN-based  implementation  of  the 
algorithm  called  BIFDE  was  designed  to  provide  a  numerical  method 
for  resolving  generic  and  1st  order  nongeneric  bifurcations  in 
functional  differential  equations  (not  necessarily  of  delay 
-difference,  or  even  finite  delay  type) .  The  code  has  been  checked 
on  specific  equations  from  the  literature,  and  has  been  shown  to  be 
useful  in  analyzing  a  relatively  unstudied  model  of  chugging  in 
liquid-propellant  fuel  rockets.  See  (1],  [2],  and  [3]. 

A  second  numerical  implementation  investigates  the  feasibility  of 
executing  the  algorithm  with  symbolic-manipulation  software.  The 
use  of  MACSYMA  for  such  purposes  has  been  undertaken,  and  results 


indicate  that  this  algorithm  succeeds  where  symbolic  implementation 
of  alternate  approaches  {eg,  center  manifold  approximation  followed 
by  the  use  of  Poincare  transforms)  have  been  unsuccessful.  The 
investigations  involve  both  the  algebraic  calculation  of  the 
bifurcation  function  for  specific  classes  of  delay-difference 
equations,  and  the  derivation  of  general  formulae  (applicable  to 
general  FDE)  for  the  resolution  of  the  bifurcation  structure  at 
points  of  second  order  nongenericity.  See  [4] . 

The  success  of  the  above  algorithm  points  out  the  need  for 
effective  means  of  computing  the  associated  bifurcation  data.  An 
important  step  in  that  direction  is  the  determination  of  the  global 
convergence  properties  of  various  rootfinding  technic[ues  when 
applied  to  quasipolynomials  of  the  type  that  arise  as 
characteristic  equations  for  linearized  problems.  Such  a  study  was 
initiated  with  a  grant  of  CPU  time  on  the  Minnesota  Supercomputer 
Center  Cray  -  2.  While  fractal  convergence  basin  boundaries  where 
anticipated,  their  precise  form  was  found  to  share  certain 
geometric  characteristics  independent  of  the  rootfinding  method  in 
use.  Comparison  of  Newton's  method  to  various  third  order  methods 
(eg.,  Halley's)  indicate  that  the  global  convergence  properties  can 
vary  greatly  from  method  to  method.  Description  of  the  results  and 
a  vector  algorithm  for  their  generation  are  given  in  [5]  and  [6]. 

The  work  involving  global  tracking  techniques  has  centered  on  the 
design  of  a  general  purpose  FORTRAN  -  based  code  for  following 
periodic  solutions  in  one  -  parameter  families  of  FDE,  identifying 
secondary  bifurcations,  and  following  certain  branches.  A  general 
purpose  code  for  the  analysis  of  scalar  delay  -  difference 
equations  was  completed  and  applied  to  various  equations  from  the 
literature.  See  [7] .  A  redesign  and  extension  of  the  code  to 
certain  systems  has  been  completed  as  well  [3] .  A  first  step  in 
designing  a  parallel  processing  version  of  the  code  was  taken  in 
[3]  . 

In  addition  to  the  FDE  taken  from  the  literature  to  test  the 
effectiveness  of  these  methods,  work  has  been  initiated  in  the 
analysis  of  a  class  of  new  problems  from  modern  electrodynamics. 

In  particular,  the  general  equations  of  motion  of  charged  particles 
(including  relativistic  and  radiation  effects)  lead  —  depending  on 
the  simplifying  assumptions  made  —  to  either  retarded  FDE,  FDE  of 
neutral  type,  or  retarded  FDE  with  state-dependent  time  delays. 

The  model's  derivation  and  the  results  of  certain  simulation 
studies  are  given  in  [9]  and  [10]. 


PROJECT  PUBLICATIONS; 

*  [1]  BIFDE:  Software  for  the  Investigation  of  the  Hopf  Bifurcation 
Problem  in  Functional  Differential  Equations  (with  N.  Aboud  and  A. 
Sathaye) ,  Proceedings  for  the  27th  IEEE  Conference  on  Decision  and 
Control,  to  appear. 

[2]  "Contributions  to  the  Computer-Aided  Analysis  of  Functional 
Differential  Equations,"  Master's  Thesis  by  N.  Aboud,  University  of 
Minnesota,  August,  1988. 

[3]  Periodicity  in  a  Model  of  Liquid-Propellant  Fuel  Rockets,  in 
preparation. 

[4]  "Symbolic  Hopf  Bifurcation  Calculations  for  Functional 
Differential  Equations,"  Master's  Thesis  by  J.  Franke,  University 
of  Minnesota,  in  progress. 
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[5]  Graphical  Analysis  of  Newton's  Method,  University  of 
Minnesota-Duluth  Mathematics  Technical  Report  #  88-4,  J. 
Reichenborn,  1988. 


[6]  Fractal  Convergence  Basins  for  Rootfinding  Techniques  Applied 
to  Transcendental  Equations,  University  of  Minnesota-Duluth 
Mathematics  Technical  Report,  E.  Swieringa,  in  progress. 

*  [7]  Computer-Aided  Analysis  of  Periodic  Solutions  of  Functional 
Differential  Equations,  Differential  Equations  and  Applications,  to 
appear . 

[8]  Parallel  Scientific  Programming  on  the  Encore  Multimax, 
University  of  Minnesota  -  Duluth  Mathematics  Technical  Report,  D. 
Kingsley,  in  progress. 

[9]  Time  Delays  in  Models  from  Electrodynamics,  (with  N.  Jahren 
and  T.  Jordan),  in  progress. 

[10]  FDESLV:  An  ODE-based  Software  Package  for  the  Simulation  of 
Functional  Differential  Equations,  University  of  Minnesota  -  Duluth 
Mathematics  Technical  Report  #88-3,  K.  Hill,  1988. 


*  6  preprints  enclosed  with  final  technical  report. 
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COMPUTER-AIDED  ANALYSIS  OF 
PERIODIC  SOLUTIONS  OF 
FUNCTIONAL  DIFFERENTIAL  EQUATIONS 

HARLAN  W.  STECH 

Department  of  Mathematics  and  Statistics 
University  of  Minnesota,  Duluth,  Minnesota  55812  USA 


I.  INTRODUCTION 

Computational  techniques  have  played  an  important  role  in  the  advancement  of 
our  understanding  of  the  qualitative  nature  of  solutions  to  functional  differential 
equations.  Simulation  studies  often  provide  compelling  evidence  for  the  asymp¬ 
totic  behavior  of  models  too  complex  for  direct  analysis.  Stable  equilibria,  periodic 
or  quasi-periodic  solutions,  or  chaotic  dynamics  are  often  more  easily  “observed" 
through  simulation  than  rigorously  proven  to  exist. 

As  important  as  simulation  studies  are,  they  are  often  restricted  in  their  effec¬ 
tiveness  by  the  fact  that  they  (in  all  but  elementary  situations)  only  identify  the 
numerical  analogue  of  “asymptotically  stable”  solutions  of  the  associated  differential 
equations.  It  is  clear  that  unstable  solutions  play  an  integral  role  in  understanding 
the  underlying  structure  of  chaotic  dynamics.  While  computational  methods  now 
exist  for  the  analysis  of  unstable  Hopf  bifurcations  in  autonomous  functional  differ¬ 
ential  equations  (see  [ij,  [3l,  [5j),  the  information  so  obtained  is  of  limited  value  in 
understanding  large  periodic  orbits. 

In  the  case  of  autonomous  ordinary  differential  equations,  the  code  AUTO  of  E. 
Doedel  [2j  assists  in  the  computation  of  ljuge  periodic  orbits  (regardless  of  their 
stability  type)  in  parameter  -  dependent  equations.  The  code  allows  (among  many 
things)  the  tracking  of  parameter  -  dependent  equilibria,  location  of  bifurcation 
points,  tracking  parameter  -  dependent  periodic  solutions  arising  through  Hopf 
bifurcations,  calculation  of  Floquet  multipliers  and  secondary  bifurcation  points, 
and  the  tracking  certain  secondary  bifurcations.  Clearly,  it  would  be  desirable  to 
have  at  one’s  disposal  a  general  purpose  code  that  can  accomplish  these  same  tasks 

AMS(MOS)  Subject  Cla«iflcation:  MCSfl,  34K99. 

This  work  wm  partially  supportad  by  grants  AFOSR  #86-0071  and  874)268,  NSF  DMS-87014S6, 
and  a  grant  from  th«  Minnasota  Supareomputar  Cantar. 


2 


’’  ?f.  STECH 


for  systems  of  autonomous  functional  differential  equations.  The  purpose  of  this 
paper  is  to  report  on  recent  progress  made  in  that  direction. 

We  discuss  a  modest  "first  step”  in  the  development  of  a  general  purpose  code 
for  the  numerical  computation  of  the  structure  of  periodic  orbits  in  parameter  - 
dependent  functional  differential  equations.  In  particular,  a  FOBTRAN  code  has 
been  recently  developed  for  the  general  scalar  delay  difference  equation 

=  /(a;^(0.  *(*  -  1)).  (1-1) 

where  a  is  a  real  (bifurcation)  parameter. 


n.  NUMERICAL  METHODS 


The  model  problem  (1.1)  has  been  chosen  as  one  of  clear  importance  in  the  math¬ 
ematical  literature,  as  well  as  one  through  which  one  can  judge  the  feasibility  of 
developing  an  analogous  code  for  general  systems  of  functional  differential  equations 
(not  necessarily  of  delay  -  difference  form).  While  (1.1)  is  of  special  structure, 
the  numerical  algorithms  employed  have  made  little  reliance  on  that  fact,  so  as  to 
provide  a  more  valid  prediction  of  the  computational  complexity  to  be  expected  for 
the  general  functional  differential  equation 


i'(t)  =/(a;xt),  (2.1) 

where  for  fixed  o,  /  :  C([— l,0j;i2'*)  -*  5'*,  the  usual  space  of  continuous  functions. 

The  development  of  a  general  purpose  code  for  the  tracking  of  periodic  solutions 
can  be  roughly  divided  into  3  numerical  issues. 


1.  Numerical  Approximation  of  Periodic  Solutions 

It  should  be  stressed  here  that  both  unstable  and  stable  periodic  solutions  are  of 
equal  importance.  Thus,  simulation  techniques,  which  would  lend  itself  to  finding 
stable  orbits,  are  not  used.  As  in  the  case  of  .411X0,  collocation  methods  have  been 
chosen. 

In  particular,  one  first  scales  time  in  (l.l)  so  that  a  periodic  solution  of  (l.l)  of 
period  T[—  1/w)  is  transformed  into  a  periodic  solution  of  period  1  in 

wx'{t)  —  /(a;i(t),  z{t  —  u>))  =  0  (2.2) 

A  partition  of  the  interval  [0.,  1.1  is  then  made  into  N  subintervals  of  equal  length. 
The  1-periodic  solution  i(t)  of  (2.2)  is  then  approximated  in  terms  of  periodic 
S-splines  of  order  k  (3  <  4  <  9).  Using  the  superscript  N  to  denote  this  approxi¬ 
mation,  we  have 

=  -  (2-3) 

where  the  coefficients  e,-,  i  =«  1,2,...  ,N  are  to  be  determined.  One  now  imposes 
N  collocation  constraints  on  by  requiring  (2.2)  hold  at  N  collocation  nodes 
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(one  taken  from  each  of  the  N  subintervais.)  Additionally,  one  must  impose  a  phase 
shift  constraint  in  order  to  factor  out  the  indeterminacy  due  the  fact  that  (2.1)  is 
autonomous.  (Thus  a  phase  shift  of  any  periodic  solution  generates  another.)  The 
resulting  syston  of  + 1  equations  in  jV  +  2  unknowns  a,  Ci,  t  =  1,2,...  liV,  and 
to  is  then  expected  to  possess  1-parameter  families  of  solutions.  Thus,  one  is  lead 
to  the  consideration  of: 

2.  One-Parameter  Curve  Tracking  Techniques 

Methods  for  computing  one-parameter  families  of  solutions  to  nonlinear  systems  are 
now  well  known.  It  is  most  convenient  to  visualize  the  family  locally  parameter¬ 
ized  in  terms  of  arc  length,  s.  Initial  points  on  the  solution  curve  can  be  obtained 
from  Hopf  bifurcation  data  near  equilibria,  or  from  spline  approximations  to  stable 
periodic  orbits  as  observed  through  simulation. 

Given  two  nearby  points  on  the  periodic  family,  one  “predicts”  another  by  linear 
extrapolation  to  a  point  at  some  prescribed  distance,  ds,  along  the  curve.  One  now 
imposes  a  final  scalar  constraint  so  as  to  insure  that  the  next  sought-after  point  on 
the  solution  curve  is  ds  arclength  units  from  the  previous  point  of  the  curve.  The 
resulting  N  +  2  equations  in  iV  -f-  2  unknowns  are  then  solved  by  Newton’s  method 
(with  numerically  determined  Jacobian). 

The  current  code  is  designed  to  implement  stepsize  selection  automatically.  Input 
data  states  the  maximum  allowable  stepsize,  as  well  as  the  optimal  and  minimum 
choices.  If  Newton’s  method  has  not  converged  after  a  certain  number  of  iterations 
(prescribed  as  input),  the  stepsize  is  reduced  and  continiution  b  attempted  again.  If 
convergence  faib  again,  the  process  b  repeated  until  either  convergence  b  obtained 
or  the  stepsize  b  reduced  to  a  value  smaller  than  the  selected  minimum  value  of  ds. 
The  latter  situation  results  in  program  termination. 

3.  Floquet  Multiplier  Approximation 

Identification  of  the  Floquet  multipliers  for  the  approximated  periodic  solution  b 
important  both  for  determining  the  solution’s  stability,  and  for  identification  of 
points  of  secondary  bifurcation.  These  multipliers  are  the  eigenvalues  of  the  lin¬ 
earized  Poincare  map  for  the  T-periodic  orbit  r(t)  of  (1.1).  It  b  convenient  to 
consider  the  linearized  equation 

2'(0  =  Dif{x{t),  x{t  -  l))z(t)  +  Dzfixit),  x{t  -  l))z{t  -  1)  (2.4) 

(suppressing  the  dependence  on  a)  on  the  phase  space  Ir*[— l,0j  x  R,  rather  than 
C([— l,0j;  iZ).  Without  loss  of  generality,  we  may  assume  T  >  1.  Thus,  the  time  T 
map  for  thb  equation  b  knowm  to  be  compact.  It  b  approximated  by  computing 
the  action  of  the  map  on  a  finite  dimensional  approximation  to  the  phase  space. 
In  particular,  one  niunerically  computes  the  solutions  of  the  initial  value  problems 
for  the  linearized  equation  with  initial  conditions  taken  to  be  one  of  the  Af  +  1 
characterbtic  functions  of  the  intervab  [— h*,  —k{i  —  1)], :  =  1, . , .  ,.Vf,  (k  =  l.JM), 
and  {0}  for  i  =  0.  Using  the  inner  product  structure  of  the  phase  space,  solutions 
on  [r  —  l.,T]  are  then  projected  onto  the  M  +  1  dimensional  subspace'  spanned  by 
the  above  initial  conditions. 
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As  a  result,  one  obtains  a  linear  map  on  the  eigenvalues  of  which  are 

expected  to  correspond  to  Flo<iuet  multipliers  for  the  periodic  solution.  Clearly, 
the  accuracy  of  the  approxiination  scheme  relies  both  on  the  accuracy  of  the  ap- 
proodmation  to  the  periodic  solution  of  (1.1),  as  well  as  the  dimension  of  finite 
dimensional  approxiznation  to  the  phase  space.  It  should  be  remarked  that  (2.1) 
being  autonf^moos  implies  that  1.  should  be  a  Floquet  multiplier,  thus  providing  a 
convenient  symptom  of  the  overall  accuracy  of  computed  multipliers. 


m.  TWO  EXAMPLES 

It 

Wright’s  Equation 

x'{t)  =  —ax{t  —  1)[1.  +  s(01.  o  >  0  (3.1) 

satisfies  the  standard  conditions  for  a  Hopf  bifurcation  at  a  =(^2.  It  is  well  known 
that  the  bifurcation  is  generic,  supercritical,  and  locally  orbitally  asymptotically 
stable  (see  [4]  and  the  references  therein.)  The  family  of  periodic  orbits  are  easily 
calculated  by  the  code  described  above.  Figure  3.1  depicts  the  bifurcation  diagram 
obtained  from  75  ctirve  tracking  steps  performed  on  the  collocation  approximation 
(5th  order  splines)  with  iV  =  28.  The  vertical  axis  is,  roughly  speaking,  the  L‘  norm 
of  the  periodic  orbit  (normalised  to  period  1.) 

Floquet  multiplier  approximation  is  reserved  (in  this  case)  to  every  10th  tracking 
step.  Because  the  multipliers  are  known  to  cluster  at  the  origin  of  the  complex 
plane,  we  tabulate  only  those  multipliers  with 
dimension  of  the  phase  space  approximation  Is 

step  multipliers  larger 


0.9104 

0.6521 


10 

0.9988 

20 

0.9991 

30 

0.9992 

40 

0.9995 

50 

1.0007 

60 

1.0085 

70 

1.0346 

modulus  larger  than  one  half.  The 
=',41. 


than  .5 


.UJx. 


Observe  that  (except  for  the  numerical  equivalent  of  the  unit  multiplier)  all  mul¬ 
tipliers  lie  within  the  unit  circle,  thus  confirming  the  stability  type  of  the  periodic 
orbit.  Note  also  that  accuracy  of  the  unit  multiplier  degenerates  as  the  periodic 
orbit  grows  in  size.  This  can  be  explained  by  the  fact  that  the  periodic  orbits  to 
Wright’s  equation  possess  regions  of  fast  variation  of  i.  Accuracy  of  the  approxi¬ 
mation  scheme  can  be  maintained  if  one  increases  the  collocation  parameter  N  as 
one  tracks  along  the  curve. 

The  above  computation  was  performed  on  the  four  -  processor  Cray-2  of  the 
University  of  Minnesota  Supercomputer  Center,  and  required  47.75  seconds  CPU 
time.  Based  on  benchmark  comparisons  with  a  VAX  11/750,  one  would  expect  the 
computation  to  take  approximately  25  minutes  on  this  more  common  machine. 
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a  final  example,  we  consider  the  delayed  negative  feedback  equation 

x'(t)  =  —2ax{t)  —  2bx{t  —  l)/[l.  +  —  1)],  (3.3) 

where  a  and  b  are  positive  parameters.  If  one  fixes  a  =  1.00,  then  the  spectral 
conditions  for  a  Hopf  bifurcation  firom  x  —  Q.  are  satisfied  at  6  =  1.5219.  Due  to 
the  choice  of  nonlinearity  the  bifurcation  is  not  generic.  However,  the  stability 
algorithm  of  [5]  can  be  applied  to  show  the  periodic  orbits^  exist  supercritically, 

ar to  be  locally  orbitally  asymptotically  stable.  - - 

Figure  3.2  depicts  the  computed  bifurcation  diagram.  Fifth  order  splines  are  used, 
with  N  and  M  varied  so  as  to  maintain  acceptable  accuracy  of  the  unit  multiplier. 
The  period  along  the  stable  leg  of  the  primary  Hopf  bifurcation  is  approximately 
2.743.  The  first  secondary  bifurcation  (at  b  =  1.953)  corresponds  to  a  simple  branch, 
with  a  Floquet  multiplier  leaving  the  unit  disk  at  1.  on  the  real  axis.  The  stability 
of  the  periodic  orbit  is  observed  to  switch  to  this  secondary  family.  This  secondary 
family  itself  looses  stability  to  a  stable  period  doubling  orbit  of  period  5.2961  at 
b  =  2.355.  These  periodic  orbits  loose  stability  at  another  period  doubling  (near 
6  =  2.453) ,  and  the  tracking  is  terminated  when  the  resulting  period  10.569  orbits 
undergo  yet  another  period  doubling  bifurcation  near  b  =  2.479  to  an  orbit  of  period 
21.113.  See  Figm  e  3.3  for  a  simultaneous  plot  of  existing  periodic  orbits  near  2.470. 
Dashed  curves  denote  unstable  orbits. 


IV.  SUMMARY 

The  examples  of  the  previous  section  indicate  the  types  of  information  that  can  be 
now  routinely  obtained  through  the  described  code.  However,  the  benchmark  data 
for  Wright’s  equation  make  it  clear  that  additional  work  must  be  done  if  a  code  of 
this  type  is  to  be  effectively  implemented  on  standau’d  mainframe  computers. 

Ongoing  research  concerns  the  effects  of  variable  mesh  collocation  methods  and 
alternate  Floquet  multiplier  approximation  methods  on  improving  the  performance 
of  the  above  code.  Tayloring  a  code  for  delay  difference  equations  is  expected  to 
improve  code  performance,  as  well. 
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